Dynamics of the Markov Time Scale of Seismic Activity May Provide 
a Short-Term Alert for Earthquakes 
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We propose a novel method for analyzing precursory seismic data before an earthquake 
that treats them as a Markov process and distinguishes the background noise from real 
fluctuations due to an earthquake. A short time (on the order of several hours) before 
an earthquake the Markov time scale Im increases sharply, hence providing an alarm 
for an impending earthquake. To distinguish a false alarm from a reliable one, we 
compute a second quantity, Ti, based on the concept of extended self-similarity of the 
data. Ti also changes strongly before an earthquake occurs. An alarm is accepted if 
both Im and Ti indicate it simultaneously. Calibrating the method with the data for one 
region provides a tool for predicting an impending earthquake within that region. Our 
analysis of the data for a large number of earthquakes indicate an essentially zero rate 
of failure for the method. 



Earthquakes are complex phenomena.^ Ahhough still 
subject to some debate, precursory anomalies, such as 
changes in the seismic recordings, and anomalous vari- 
ations in the chemical, hydrological, and electromag- 
netic properties of the area in which earthquakes oc- 
cur, usually precede a large earthquake. One school 
of thought believes that the anomalies occur within days 
to weeks before the main shock, but probably not much 
earlier,^ and that the spatial precursory patterns develop 
at short distances from impending large earthquakes. A 
second school believes that the anomalies may occur up 
to decades before large earthquakes, at distances much 
larger than the length of the main shock rupture, a con- 
cept closely linked to the theory of critical phenomena^^^ 
which was advocated^'^'^ as early as 1964 with a report 
documenting the existence of long-range correlations in 
the precursors. Knopoff et al.^ reported recently the exis- 
tence of long-range spatial correlations in the increase of 
medium-range magnitude seismicity prior to large earth- 
quakes in California. 

Pursuing a model of rock rupture and its relation 
with critical phenomena and earthquakes,^ a method 
of analysis was introduced^'^ that, for certain values 
of its parameters, led to a power law (typical of criti- 
cal phenomena) for the system's time-to-failure. Several 
groups^" proposed percolation^^ and hierarchical models 
of damage/rupture prior to an earthquake. In particular, 
Sahimi et al.^^ proposed a connection between percola- 



tion, the spatial distribution of earthquakes' hypocen- 
ters, and rock's fracture/fault networks. Sornette and 
Sammis^^ developed a theory according to which the 
power law that describes the accelerated seismicity close 
to a large earthquake is accompanied by log-periodic cor- 
rection terms, which were shown^'^ to also exist in the 
power law that describes the increase in the energy that 
rock releases as it undergoes fracturing. Such ideas were 
further developed by Huang et al.,^^ with empirical ev- 
idence provided by Bowman et al.,^'^ and view a large 
earthquake as a temporal singularity in the seismic time 
series, resulting from the collective behavior and accu- 
mulation of many previous smaller-size events. In this 
picture, as the stress on rock increases, micro-ruptures 
develop that redistribute the stress and generate fluctu- 
ations in it. As damage accumulates, the fluctuations 
become spatially and temporally correlated, resulting in 
a larger number of significantly-stressed large domains. 
The correlations accelerate the spatial smoothing of the 
fluctuations, culminating in a rupture with a size on the 
order of the system's size, and representing its final state 
in which earthquakes occur. Numerical^^ and empirical^" 
evidence for this picture indicates that, similar to criti- 
cal phenomena, the correlation length of the stress-field 
fluctuations increases significantly before a large earth- 
quake. Notwithstanding the evidence, proving or refut- 
ing the notion of earthquakes as a critical phenomenon 
entails replacing the proxies, used for checking the build- 



up of the coopcrativity that leads to large earthquakes, 
by a direct measure of the dynamic evolution of the stress 
field. Unfortunately, such a procedure is far beyond the 
present technical abilities. 

A theory of earthquakes should predict, (1) when and 
(2) where they occur in a wide enough region. It should 
also be able to (3) distinguish a false alarm from a re- 
liable one. In this paper, we propose a method for pre- 
dicting earthquakes which possesses the three features. 
The method estimates the Markov time scale (MTS) tj^j 
of a seismic time series - the time over which the data 
can be represented by a Markov process. ^^"^^ As the 
seismic data evolve with the time, so also does tnj- We 
show that the time evolutioon of Im provides an effec- 
tive alarm a short time before earthquakes. The method 
distinguishes abnormal variations of ti^i before the arrival 
of the P- waves, hence providing enough of a warning for 
triggering a damage/death-avoiding response prior to the 
arrival of the more damaging S- waves. 

The method first checks whether the seismic data 
follow a Markov chain and, if so, measures the func- 
tion MTS Characterization of the statisti- 
cal properties of fluctuations of n measured quanti- 
ties of the stochastic process x{t) requires evaluation 
of the joint probability distribution function (PDF) 
Pn{xi,ti; ■ ■ ■ ; Xn,tn)- If the data are a Markov pro- 
cess, then, Pn = U"~^p{xi+i,ti+i\xi,ti)p{xi,ti), where 
p{xi+i,ti^i\xi,ti) are conditional probabilities such that 
the Chapman-Kolmogorov (CK) equation, 

p{x2,t2\xi,ti) =^ J dx3 p{x2,t2\x3,t3)p{x3,t3\xi,ti) (1) 

holds for any ^3 in ti < < 12- The validity of 
the CK equation for describing the process is checked 
by comparing the directly-evaluated p(x2,t2\xi,ti) with 
the those calculated according to right side of Eq. 
(1). To determine Im for the data we compute for 
given xi and X2 the quantity, Q = \p{x2,t2\xi,ti) — 
J dx3p{x2, t2\x3, <3)p(a;3, tsjcci, ti)|, in terms of, for exam- 
ple, ^3 — <i. In practice, we take ti = and ^3 = ^t2, 
and vary t2. Plotting Q versus t3 produces the position 
of tM in the limit Q 0.^^ 

Our analysis of seismic data (see below) indicates that 
the average Im for the uncorrected background seismic 
time series is much smaller than that for earthquakes data 
(P-wave plus S-wave). Thus, at a certain time before an 
earthquake, Im rises significantly and provides an alarm 
for the earthquake. As we show below, the alert time ta is 
on the order of hours, and depends on the earthquake's 
magnitude M and the epicenter's distance d from the 
data-collecting station(s). 

The sharp rise in tM at the moment of alarm is, in some 
sense, similar to the increase in the correlation length 
^ of the stress-field fluctuations in the critical phenom- 
ena theories of earthquake, since tM is also the time over 
which the events leading to an earthquake are correlated. 



Therefore, just as the correlation length ^ increases as the 
catastrophic rupture develops, so also docs tM- However, 
whereas it is exceedingly difficult to directly measure ^, 
tjv/ is computed rather readily. Moreover, whereas ^ is 
defined for the entire rupturing system over long times, 
tjvf is computed online (in real time), hence reflecting the 
correlations of the most recent events that are presum- 
ably most relevant to an impending earthquake. 

To distinguish a false alarm that might be indicated 
by tM from a true one, we use a second time-dependent 
function that we compute based on the extended self- 
similarity (ESS) of the seismic time series. ^"''■^'^ The ESS 
is characterized by Sp, a structure function of order p, 
defined by 

Sp{r) = {\xit + r) - x{tW) - Mt + r) ~ x{t)f)<^ (2) 

where r is the lag (in units of data points). The first 
nontrivial moment (beyond the average and variance) of 
a distribution is S3, and because for a Gaussian process, 

^ ^P' ^^^'■^ deviations from this relation represent non- 
Gaussian behavior. It is also well-known^^^^^ that the 
moments Sp with p < 1 contain information on frequent 
events in a time series. Prior to an earthquake the num- 
ber of frequent events (development of cracks that join 
up) suddenly rises, indicated by a sudden change in Sp 
with p < 1. We observe that the starting point of Sp{T) 
{p < 1) versus 5*3 (t) is different for different type of data 
set.^^'^^ To determine the distance form the origin we de- 
fine the function Ti = T{t = 1) ^ [Sq i{t = 1) + S'|(t = 
l)]i/2^ Close to an earthquake the function Ti{t), also es- 
timated online, suddenly changes and provides a second 
alert. Its utility is due to the fact that it is estimated 
very accurately even with very few data points, say 50, 
hence enabling online analysis of the data collected over 
intervals of about 1 second. Thus, even with few data 
points, the method can detect the change of correlations 
in the incoming data. For example, for correlated syn- 
thetic data with a spectral density 1//^"^^, one obtains 
Ti = -7a + 7. 

We have analyzed the data for vertical ground 
velocity Vz(t) for 173 earthquakes with magnitudes 
3.2 < M < 6.3 that occurred in Iran between 
28°N and 40°N latitude, and 47°E and 62.5°E 
longitude, between January 3 and July 26, 2004. 
Recorded by 14 stations, the data can be accessed at 
http://www.iiees.ac.ir/bank/bank_2004.html. The fre- 
quency was 40 Hz for 2 of the stations and 50 Hz for 
the rest. The vertical ground velocity data were ana- 
lyzed because with our method they provide relatively 
long (on the order of several hours), and hence useful, 
alarms for the impending earthquakes. Fourty (discrete) 
data points/second are recorded in the broad-band seis- 
mogram for the vertical ground velocity x{t) = Vz- To 
analyze such data and provide alarms for the area for 
which the data are analyzed, we proceed as follows. 




FIG. 1. Time-dependence of Ti and tM for a recent earth- 
quake of magnitude 6.3 in northern Iran, and their comparison 
with the vertical ground velocity data Vz{t). tM is in number 
of data points (the frequency at the station is 40 Hz), Ti is 
dimensionless, while 14 (t) is in "counts" which, when multi- 
plied by a factor 1.1382 x 10""^, is converted to /xm/sec. The 
sensors were (broad-band) Guralp CMG-3T that collect data 
in the east- west, north-south and vertical directions. The 
thresholds are Imc = 5.6 and Tic = 0.88. 

(1) The data are analyzed in order to check whether 
they follow a Markov chain [the directly-computed 
p{x2,t2\xi,ti) must be equal to the right side of Eq. (1)]. 
(2) The MTS tmit) of the data are estimated by calcu- 
lating the residual Q of the CK equation (see above). 
For long-enough data series (10'^ data points or more) 
the function tmit) are estimated where Q 0, but for 
shorter series the minimum in Q provides estimates of 
tM{t). (3) Ti(<) is computed for the same data. To com- 
pute S'p(t) (we used p = 1/10) the data x{t) are nor- 
malized by their standard deviation, hence making Ti 
dimensionless. (4) Steps (l)-(3) are repeated for a large 
number of previously-occurred earthquakes of size M at 
a distance d from the station, referred to as (M, d) earth- 
quakes. Earthquakes with M < Mc and d > dc are of no 
practical importance and are ignored (we used Mc = 4.5 
and dc = 150 km). (5) Define the thresholds tnjc and 
Tic such that for tM > thic and Ti > Tic one has an 
alert for an earthquake [M > Mc,d < dc). If <a/c and 
Tic are too large no alert is obtained, whereas one may 
receive useless alerts if they are too small. By comparing 
the data for all the earthquakes with M > Mc regis- 
tered in a given station, tMc and Tic for the earthquakes 
are estimated. (6) Real-time data analysis is performed 
to compute the function t]\j{t) and Ti{t). An alarm is 
turned on if > ^Mc and Ti > Tic simultaneously. 
When the alarm is turned on, it indicates that an earth- 
quake of magnitude M > Ale at a distance d < dc is 



going to occur. The procedure can be carried out for any 
station. The larger the amount of data, the more precise 
the alarm will be. 

Figure 1 presents Ti{t) and tmit) for an M ~ 6.3 
earthquake, occurred on May 28, 2004 at 12:36 am in 
Baladeh at (36.37N, 51.64E, depth 28) in northern Iran. 
The data were collected at Karaj station (near Tehran, 
Iran) at a distance of 74 km from the epicenter, and a 
depth of 70 m. The earthquake catalogue in the inter- 
net address given above indicates that, for several days 
before the main event, there was no foreshock in that re- 
gion. Thus, Ti and tM provided a seven hour alarm for 
the Baladeh earthquake. Since the data used for comput- 
ing tM and Ti were, respectively, in strings of 200 and 50 
points, there is no effect of the events before they were 
collected and, hence, the patterns in Fig. 1 reflect the 
events taking place in the time period in which the data 
were collected. 

To estimate the alert times ta, which are on the or- 
der of hours, we carried out an analysis of online data 
for 14 stations in Iran's broad-band network (the sensors 
are Guralp CMG-3T broad-band), analyzing the verti- 
cal ground velocity data. Our analysis indicates that ta 
depends on M, being small for low M, but quite large 
for large M. Using extensive data for the Iranian earth- 
quakes with M > 4.5 and d < 150 km, we have ob- 
tained an approximate relation for the broad-band sta- 
tions, shown in Figure 2 and represented by 

logf„ = -1.35-^2.4 log M , (3) 

where ta is in hours. The numerical coefficients of Eq. 
(3) for each area should be estimated from the data col- 
lected for that area. The above analysis can clearly be 
extended to all the stations around the world. This is 
currently underway for Iran's network. For an earth- 
quake of magnitude M = 4.5, Eq. (3) predicts an alert 
time of about 2 hours. Thus, if, for example, three hours 
after the alarm is turned on, the earthquake has not still 
happened, we know that the magnitude of the coming 
earthquake is M > 5.7. 

In summary, we have proposed a new method for ana- 
lyzing seismic data and making predictions for when an 
earthquake may occur with a magnitude M > Mc at a 
distance d < dc- The method is based on computing 
the Markov time scale and a quantity Ti calculated 
based on the concept of extended self-similarity of the 
data, and monitoring them online as they evolve with the 
time. If the two quantities exceed their respective critical 
thresholds tcM and Td, estimated based on analyzing the 
data for the previously-occurred earthquakes, an alarm 
is tuned on. We are currently utilizing this method for 
Iran's stations. To do so, we calibrate the method with 
the data for the stations in one region (i.e., estimate tcM 
and Tci for distances d < dc). If in a given region there 
is a single station, then once the online-computed tM 



and Ti exceed their eritical values, the alarm is turned 
on. If there are several stations, then once they declare 
that their tm and Ti have exceeded their thresholds, the 
alarm is turned on. If after about 2 hours, no earthquake 
has occurred yet, then we know that the magnitude of 
the incoming earthquake will be greater Mc = 4.5 at a 
distance d < dc- 

Over the past two years, the method has been utilized 
in the Iranian stations. Our analysis indicates that the 
method's failure rate decreases to essentially zero when 
Im and Ti provide simultaneous alarms. That is, practi- 
cally every earthquake that we have considered, including 
those that have been occurring while we have been per- 
forming online analysis of their incoming data and pro- 
viding alarms for them (with M > Mc), was preceded by 
an alarm. Of all the earthquakes that we have analyzed 
so far, the method has failed in only two cases. In our 
experience, if after 10 hours no earthquake occurs, we 
count that as a failed case. However, as mentioned, we 
have so far had only two of such cases. 

Finally, it must be pointed out that the most accurate 
alarms are obtained from stations that receive data from 
depths of > 50 m, and are perpendicular to the active 
faults that cause the earthquake, since they receive much 
more correlated data for the development of the cracks 
than any other station. 
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